Signatures of Dynamical Heterogeneity in the Structure of Glassy Free-energy Minima 
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From numerical minimization of a model free energy functional for a system of hard spheres, we 
show that the width of the local peaks of the time-averaged density field at a glassy free-energy 
minimum exhibits large spatial variation, similar to that of the "local Debye-Waller factor" in 
simulations of dynamical heterogeneity. Molecular dynamics simulations starting from a particle 
configuration generated from the density distribution at a glassy free-energy minimum show similar 
spatial heterogeneity in the degree of localization, implying a direct connection between dynamical 
heterogeneity and the structure of glassy free energy minima. 



PACS numbers: 

Observation of dynamical heterogeneity, both in ex- 
periments [l[ and in simulations [2], has been a signif- 
icant step in the attempt to understand the behaviour 
of glass forming liquids. The degree of spatial variation 
of the "propensity for motion" 0] of individual parti- 
cles, defined as the mean-square displacement of a par- 
ticle from its initial position, provides a direct measure 
of dynamical heterogeneity. The local Debye-Waller fac- 
tor (short-time rms displacement from the average po- 
sition) of individual particles in the initial configuration 
has been found [ij to be strongly correlated with the spa- 
tially heterogeneous propensity for motion over longer 
time scales. However, subsequent work |5| has shown 
that this correlation exists only for time scales shorter 
than the a-relaxation time. Spatial heterogeneity of the 
local Debye-Waller factor has also been observed [6[ in 
simulations below the glass transition temperature. 

The physical origin of dynamical heterogeneity is not 
well-understood at present. In particular, it is not clear 
whether the occurrence of spatially heterogeneous dy- 
namics can be explained within the "free-energy land- 
scape" description [3, @, Q of glassy behavior in which 
the complex dynamics is attributed to the presence of 
a large number of "glassy" local minima of the free en- 
ergy. Density functional theory (DFT) [l(|, in which 
the free energy is expressed as a functional of the time- 
averaged local density, provides a convenient framework 
for exploring the free-energy landscape. In this descrip- 
tion, a glassy free-energy minimum is a local minimum of 
the free energy functional with a strongly inhomogeneous 
but non-periodic density distribution. The local peaks of 
the density distribution represent the time-averaged po- 
sitions of the particles, and the width of a local peak is 
analogous to the local Debye-Waller factor measured in 
simulations. In this description, the a-relaxation time 
corresponds to the time scale of transitions between dif- 
ferent glassy minima Therefore, the density dis- 
tribution at a typical glassy free-energy minimum should 
correspond to an average of the local density over a time 
scale shorter than the a-relaxation time. If this descrip- 



tion is valid, then the spatial variation of the propensi- 
ties for motion observed in simulations over time scales 
shorter than the a-relaxation time (which, as discussed 
above, is strongly correlated with the spatial variation 
of the local Debye-Waller factor) should be manifested 
in the structure of glassy free-energy minima as a simi- 
lar spatial variation of the widths of the local peaks of 
the density distribution. The observation of such spa- 
tial variation of the peak width would, therefore, pro- 
vide an explanation of dynamical heterogeneity within 
the "free-energy landscape" description. This would be 
an alternative to the recently proposed "mode-coupling" 
description [11] of dynamical heterogeneity. 

In this Letter, we have used numerical minimization 
of the Ramakrishnan-Yussouff (RY) free energy func- 
tional [10|| for a hard sphere system to study the density 
distribution at glassy free-energy minima. Using a Gaus- 
sian superposition approximation [liL [l3L [lit ] , as well as 
unconstrained numerical minimization of a discretized 
version of the free-energy functional Q, we show that 
the width of the local density peaks at glassy minima ex- 
hibits large spatial variation similar to the spatial hetero- 
geneity of the local Debye-Waller factor found in simula- 
tions 0, 0] . We have also performed molecular dynamics 
(MD) simulations starting from a particle configuration 
corresponding to a glassy minimum. The simulation re- 
sults for the rms displacement of the particles from their 
average positions are found to be strongly correlated with 
the corresponding widths of the local density peaks ob- 
tained in the DFT calculation. These results establish a 
direct connection between dynamical heterogeneity and 
the structure of glassy free-ener gy m inima. 

The RY free energy functional [l0( for a system of hard 
spheres has the form 



PF = J dv{p(T)ln(p(T)/p )-Sp(v)} 

- \JdvJ dv'C{\v-v'\)8p{r)§p{v'), (1) 
where (3 = l/(fcsT), T is the temperature, 5p(v) = p(r) — 
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po is the deviation of the time-averaged local density p(r) 
from the density po of the uniform liquid, and C(r) is the 
direct pair correlation function of a uniform hard-sphere 
liquid at density p . In Eq.JT]), we have taken the zero 
of the free energy at its uniform liquid value. In earlier 
studies [H, [H, [3] of glassy minima of this free energy, 
the local density p(r) was approximated as a superposi- 
tion of Gaussians, p(r) = Yj%=i exp[— a(r — Ri) 2 ], 
with the centers of the Gaussians, {Ri}, forming a fixed 
amorphous structure which was taken to be either a ran- 
dom close packing of hard spheres 0, [r| , or a particle 
configuration from MD simulations [14| . The free energy 
was then minimized with respect to the single variational 
parameter a. In all these calculations, glassy states with 
free energy lower than that of the uniform liquid were 
found at high densities. 

Since we want to explore the possibility of the degree 
of localization (measured by the width parameter a) be- 
ing different for different particles, we have assumed the 
following form for the density profile: 
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where a.i and r\{ are, respectively, parameters that char- 
acterize the width and the height of the peak of the den- 
sity profile at the point R^. In our calculations, initial 
values of {Ri} are taken from particle configurations ob- 
tained from MD simulations of N = 500 hard spheres. 
We then minimize the RY functional with respect to the 
2500 parameters, {Ri}, {aj}, {?7i}, and also with respect 
to the volume V at a reference liquid packing fraction 
(pi = 7T/9oc 3 /6, a being the hard sphere diameter. We use 
the Hendersen-Grundke (HG) expressio n [15 j for C(r), as 
has been done in previous calculations |l2l. Il3l. Hij . The 
minimization leads to structures similar to glassy states 
observed in simulations and experiments. In Fig[TJ the 
pair distribution function, g(r), of the set of coordinates 
{Ri} for one such minimum has been plotted - the split 
second peak of the g(r) clearly indicates that the corre- 
sponding structure is amorphous (l6| . The inset of Figfl] 
shows the dependence of the free energy of this minimum 
on the packing fraction <pi . The free-energy of the glassy 
minimum becomes lower than that of the uniform liq- 
uid (i.e. /3F becomes negative) at (pi = 0.526, which, as 
expected, is lower than the value obtained in Ref. [rJ 03] • 
In order to measure the widths of the local peaks of 
the density field at a glassy minimum, we have evaluated 
the quantities (Sr)i = (j v dr\r — Ri| 2 /j(r))2 ; where Vi 
is a small volume centered at the point R^ (position of 
the ith local peak), such that J drp(r) = 1. We find 
that, indeed, the values of {8r)i are widely distributed, 
as can be seen from Fig|2j where we have plotted the 
distribution of (Sr), scaled by its spatial average (<5r) a v, 
for a glassy minimum at <pi = 0.55. The distribution in 
Fig. [5] is qualitatively similar to those of the dynamical 
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FIG. 1: The pair distribution function g(r) of a glassy mini- 
mum at (pi = 0.526 obtained using the Gaussian superposition 
approximation and the HG form of C(r). The inset shows the 
variation of the dimensionless free energy (/3F) of the glassy 
minimum with the reference liquid packing fraction pi. 



propensity and the local Debye- Waller factor obtained in 
simulations d, 0, @]. 
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FIG. 2: Distribution of (Sr), scaled by its average value, for 
a glassy minimum at pi — 0.55 obtained using the Gaussian 
approximation and the HG C(r). The distribution of the rms 

displacement (Sr)^^ obtained from MD simulations is also 
shown for comparison. Inset: correlation between (Sr)^^ 
and (Sr) (the straight line is a guide to the eye). 



To make a direct comparison between these DFT re- 
sults and the dynamics of the system, we have performed 
MD simulations of the dynamics of the system near a 
free energy minimum. In these simulations, the values of 
{Ri} obtained in the free-energy minimization are taken 
to be the initial particle positions. A few occurrences of 
neighboring local density peaks separated by less than 
a are removed using a local relaxation procedure 



Starting from the same initial configuration but using 
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different sets of initial velocities randomly assigned from 
the Maxwell-Boltzmann distribution, we have simulated 
the dynamics for 10 5 collisions per particle, which corre- 
sponds to the middle of the plateau of the mean-square 
displacement of the particles from their initial positions. 
Using 100 such different trajectories, we have calculated 
the rms displacements, {(Sr)j }, of all the particles 
from their average positions. As shown in Fig. [2J the 
distribution of (Sr) scaled by its average over all the 
particles is nearly identical to the distribution of (Sr) ob- 
tained in the DFT calculation. 

To test whether particles that have (small) large val- 
ues of (Sr) in the DFT calculation show (small) large 
values of the rms displacement in the MD simulation, we 
have performed free-energy minimizations starting from 
15 particle configurations obtained in the MD simulation. 
To account for small differences between the initial parti- 
cle positions and the positions of the local density peaks 
at the corresponding free-energy minimum, we have di- 
vided the system into 125 "cells" of equal volume and 
averaged both the DFT and MD values of (Sr) over the 
particles lying in each cell. The inset of Fig. [^illustrates 
the correlation between these spatially "coarse-grained" 
values of (Sr) obtained from DFT and MD calculations. 
The degree of correlation is quite large, with a correla- 
tion coefficient of about 0.76. The correlation coefficient 
is larger (~ 0.85) for the subset consisting of 10% of data 
points with the highest values of (Sr) / (Sr) av ■ 

Next, we have removed the Gaussian constraint on the 
density field and obtained glassy free-energy minima us- 
ing a numerical scheme [8j for minimizing a discretized 
version of the RY functional for hard spheres. To dis- 
cretize the RY functional, a simple cubic computational 
mesh of size L 3 and with periodic boundary conditions 
is introduced. On the sites of this mesh, we define den- 
sity variables pi = p(ri)h 3 , where pfa) is the density at 
site i and h the spacing of the computational mesh. The 
free energy of Eq. ([1]) is then numerically minimized as a 
function of the density variables [8J]. 

In the numerical minimization, we use the Percus- 
Yevick (PY) expression for C(r) [18| because it leads to 
(see below) a very accurate value of the crystallization 
density. The input density field {pi} for the minimization 
is obtained using configurations from a MD simulation, 
for a mesh-spacing h ~ 0.05<r. From the minimization of 
the free energy as a function of 4.096 x 10 6 density vari- 
ables, we have again obtained local minima with glassy 
{pi}. The structure of a local minimum can be charac- 
terized by the two-point correlation function G(r) of the 
time-averaged local density variables pi at the minimum 
(this function is different from the pair distribution func- 
tion g(r) of Fig. Q]). In Fig[3l we have plotted the G(r) 
for a glassy free-energy minimum at <pi — 0.586. The 
glassy nature of the density distribution is indicated by 
the split second peak of G(r) and the positions of the two 
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FIG. 3: The correlation function G(r) of the time-averaged 
local density for a glassy minimum at (f>i = 0.586, obtained 
from unconstrained minimization using the PY C(r). The 
inset shows the variation of the free energy of the glassy min- 
imum and that of the fee crystal with the packing fraction 



sub-peaks The inset of Fig[3] shows the variation of 
the free energy with the packing fraction 0/, which indi- 
cates that the free energy of the glassy minimum crosses 
that of the liquid at <fii = 0.554 which is very close to 
the packing fraction at the ideal glass transition of mode- 
coupling theory. Also shown in the plot is the free energy 
of the fee crystalline minimum, which becomes negative 
at (pi ~ 0.497, a value very close to the known packing 
fraction at crystallization. 




FIG. 4: Distribution of (Sr) at a glassy minimum obtained 
using unconstrained minimization with the PY C(r) for (f>i — 
0.613. The distribution of the same quantity obtained using 
the Gaussian approximation with the same C(r) and the dis- 
tribution of (Sr)^^ obtained from MD simulations are also 
shown for comparison. 

Using the values of {pi} at a glassy minimum, we have 
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identified the local density peaks and calculated their 
widths {{Sr)i}, as in the calculation using the Gaus- 
sian approximation. It is clear from FigJH which shows 
the distribution of (Sr) for a minimum at <pi = 0.613, 
that glassy minima obtained using this scheme also ex- 
hibit considerable spatial heterogeneity in the values of 
(Sr). To relate these results to the dynamics, we have 
performed MD simulations, as described above, start- 
ing from the particle configuration obtained from a free- 
energy minimum found in the unconstrained minimiza- 
tion. These runs consist of 50 different trajectories, each 
for 10 5 collisions per particle. As shown in Fig. [4] the dis- 
tribution of the rms displacements obtained from these 
MD runs is quite similar to the distribution of (Sr) ob- 
tained in the DFT calculation. The means and standard 
deviations of the two distributions, 0.35cr and 0.08cr, re- 
spectively, from MD and 0.34cr and 0.09cr from DFT, are 
nearly the same. Typical values of (5r)i obtained using 
the HG C(r) are about three times smaller than those 
shown in Fig. |4] Similar unrealistically small values of 
the peak widths were obtained in earlier DFT calcula- 
tions [LJ, [lj] using the HG C(r). We expect the results 



obtained using the PY C(r) to be more reliable. 
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FIG. 5: Correlation of the (5r) for a glassy minimum at 
4>i = 0.613 (obtained using the Gaussian approximation with 

the PY C(r), see text) with the rms displacement (Sr)^^ 
obtained from MD simulations. The straight line is a guide 
to the eye. 

Since the local density peaks at the free-energy mini- 
mum obtained using the PY C (r) are much broader, the 
process of elimination of neighboring peaks separated by 
less than a causes substantial rearrangement of the parti- 
cles. For this reason, a particle-by-particle comparison of 
these DFT results with the corresponding MD results is 
not possible. To make such a comparison, we have done 
a new calculation using the Gaussian approximation, but 
this time with the PY C(r), to locate a free-energy min- 
imum in which the local density peaks coincide with the 
particle positions in the initial configuration used in the 



MD runs. In this calculation, the density profile is as- 
sumed to have the form of Eq.([2]) and the free energy 
is minimized with respect to the parameters {ai} and 
{f]i}, keeping {Ri} fixed. Fig. [5] shows the correlation 
between the {(Sr)i} obtained from this DFT calculation 

and the {(Sr)^ } obtained from MD simulations (both 
spatially "coarse-grained" using 125 cells as before). This 
plot shows a high degree of correlation, with a correla- 
tion coefficient of 0.75 (it increases to 0.87 for 10% of data 
points with the highest values of (Sr)). Also, as shown in 
FigHJ the distribution of (Sr) obtained from this Gaus- 
sian DFT calculation is similar to those obtained from 
the other calculations. These results and those shown 
in Fig. [2] establish a direct connection between dynami- 
cal heterogeneity and the structure of glassy free-energy 
minima by showing that if a small region of a free-energy 
minimum contains local density peaks with small (large) 
widths, then the particles in the same region are very 
likely to exhibit small (large) values of rms displacement 
for dynamics near this minimum. 

To summarize, using DFT, we have shown that at a 
glassy minimum of the free-energy of the hard-sphere sys- 
tem, the degree of localization of the local particle den- 
sity exhibits large spatial variation. Also, by carrying 
out MD simulations starting from a particle configura- 
tion corresponding to the density distribution at a glassy 
minimum, we have established that the time-averaged 
density distribution at the minimum is closely related to 
the dynamical heterogeneity exhibited by the system. 
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